function diff = rho_fun( gamma, region, param)

% fundamental quadratic

    if strcmp(region, 'N')
        diff = -1/2.*param.sigma^2.*gamma.^2 ...
            -(param.mu - 1/2.*param.sigma^2 + (1-param.alpha).*(param.s.*param.lambda+param.zeta)).*gamma ...
            + param.r + param.s.*param.lambda + param.zeta;
    elseif strcmp(region, 'R')
        diff = -1/2.*param.sigma^2.*gamma.^2 ...
            -(param.mu - 1/2.*param.sigma^2).*gamma ...
            + param.r;
    end

end